function rgb_img = lab2rgb(lab_img)

%Step 1
xyz_img = lab2xyz(lab_img);
lrgb_img = xyz2lrgb(xyz_img);
nrgb_img = lrgb2nrgb(lrgb_img);
rgb_img = denormalized(nrgb_img);
end
%--------------------------------------------------------------------------
function rgb_img = denormalized(nrgb_img)
r_channel = nrgb_img(:,:,1);
g_channel = nrgb_img(:,:,2);
b_channel = nrgb_img(:,:,3);

if max(max(r_channel)) < 1.0 || max(max(g_channel)) < 1.0 || max(max(b_channel)) < 1.0
  r_channel = floor(double(r_channel) * 255);
  g_channel = floor(double(g_channel) * 255);
  b_channel = floor(double(b_channel) * 255);
end

rgb_img(:,:,1) = r_channel;
rgb_img(:,:,2) = g_channel;
rgb_img(:,:,3) = b_channel;
end
%--------------------------------------------------------------------------
function nrgb_img = lrgb2nrgb(lrgb_img)
[rows cols channels]=size(lrgb_img);
T = 0.0031308;

lr = lrgb_img(:,:,1);
lg = lrgb_img(:,:,2);
lb = lrgb_img(:,:,3);

nr(lr <= T) = lr(lr <= T) * 12.92;
nr(lr > T) = (1 + 0.055) * lr(lr > T) .^ (1/2.4) - 0.055;

ng(lg <= T) = lg(lg <= T) * 12.92;
ng(lg > T) = (1 + 0.055) * lg(lg > T) .^ (1/2.4) - 0.055;

nb(lb <= T) = lb(lb <= T) * 12.92;
nb(lb > T) = (1 + 0.055) * lb(lb > T) .^ (1/2.4) - 0.055;

nrgb_img(:,:,1) = reshape(nr, rows, cols);
nrgb_img(:,:,2) = reshape(ng, rows, cols);
nrgb_img(:,:,3) = reshape(nb, rows, cols);
end
%--------------------------------------------------------------------------
function srgb_img = xyz2lrgb(xyz_img)
[rows cols channel] = size(xyz_img);
X = xyz_img(:,:,1);
Y = xyz_img(:,:,2);
Z = xyz_img(:,:,3);

X = reshape(X, 1, []);
Y = reshape(Y, 1, []);
Z = reshape(Z, 1, []);

M = [3.2406 -1.5372 -0.4986; -0.9689 1.8758 0.0415; 0.0557 -0.2040 1.0570];

temp_srgb_img = M * [X; Y; Z];
sR = temp_srgb_img(1,:);
sG = temp_srgb_img(2,:);
sB = temp_srgb_img(3,:);
sR = reshape(sR, rows, cols);
sG = reshape(sG, rows, cols);
sB = reshape(sB, rows, cols);

srgb_img(:,:,1) = sR;
srgb_img(:,:,2) = sG;
srgb_img(:,:,3) = sB;
end
%--------------------------------------------------------------------------
function xyz_img = lab2xyz(lab_img)
[rows cols channels] = size(lab_img);
l_channel = lab_img(:,:,1);
a_channel = lab_img(:,:,2);
b_channel = lab_img(:,:,3);

Xn = 0.950;
Yn = 1;
Zn = 1.089;

la = (l_channel + 16)/116 + a_channel/500;
l = (l_channel + 16)/116;
lb = (l_channel + 16)/116 - b_channel/200;

T = 0.008856^1/3;

if_la (la > T) = la(la > T) .^ 3;
if_la (la < T) = (la(la < T) - 16/116)/7.787;

if_l (l > T) = l(l > T) .^ 3;
if_l (l < T) = (l(l < T) - 16/116)/7.787;

if_lb (lb > T) = lb(lb > T) .^ 3;
if_lb (lb < T) = (lb(lb < T) - 16/116)/7.787;

x_channel = Xn * if_la;
y_channel = Yn * if_l;
z_channel = Zn * if_lb;

xyz_img(:,:,1) = reshape(x_channel, rows, cols);
xyz_img(:,:,2) = reshape(y_channel, rows, cols);
xyz_img(:,:,3) = reshape(z_channel, rows, cols);
end